%% simulate linearized model

y_ss  = oo_.mean(1);
i_ss  = oo_.mean(2);
xi_ss = oo_.mean(3);
pr_ss = oo_.mean(4);
psi_p = oo_.mean(5);
n0_p  = oo_.mean(6);
gam_p = oo_.mean(7);
k_ss  = oo_.mean(8);       %state
n_ss  = oo_.mean(9);       %state
l_ss  = oo_.mean(10);      %state
rb_ss = oo_.mean(11);      %state
rkhats_ss = oo_.mean(12);  %state
a_ss  = oo_.mean(13);      %state
c_ss  = oo_.mean(14);
h_ss  = oo_.mean(15);

%% simulate

T = 10000;                    % simulation periods
rng(1)                        % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock

% policy functions
c_sim      = zeros(1,T+1);
y_sim      = zeros(1,T+1);
i_sim      = zeros(1,T+1);
k_sim      = zeros(1,T+1);
h_sim      = zeros(1,T+1);
n_sim      = zeros(1,T+1);
l_sim      = zeros(1,T+1);
rb_sim     = zeros(1,T+1);
rkhats_sim = zeros(1,T+1);
a_sim      = zeros(1,T+1);
pr_sim     = zeros(1,T+1);
xi_sim     = zeros(1,T+1);
pib_sim    = zeros(1,T+1);
n_sim_pib  = zeros(1,T+1);

lambda_sim     = zeros(1,T+1);

ss_value = oo_.dr.ys;
dr = oo_.dr.ghx;

% set initial state
k_sim(1) = k_ss;
n_sim(1) = n_ss;
l_sim(1) = l_ss;
rb_sim(1) = rb_ss;
rkhats_sim(1) = rkhats_ss;
a_sim(1) = a_ss;

minus_profit = 0;

for t=2:T
    
    c_sim(t) = c_ss + dr(14,1)*(k_sim(t-1) - k_ss) + dr(14,2)*(n_sim(t-1)     -     n_ss) + dr(14,3)*(l_sim(t-1)-l_ss) ...
                    + dr(14,4)*(rb_sim(t-1)-rb_ss) + dr(14,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(14,6)*(a_sim(t-1)-a_ss);
    h_sim(t) = h_ss + dr(15,1)*(k_sim(t-1) - k_ss) + dr(15,2)*(n_sim(t-1)     -     n_ss) + dr(15,3)*(l_sim(t-1)-l_ss) ...
                    + dr(15,4)*(rb_sim(t-1)-rb_ss) + dr(15,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(15,6)*(a_sim(t-1)-a_ss);
    lambda_sim(t) = 1 / ( c_sim(t) - psi_p*h_sim(t)^(1+1/nu_p) / (1+1/nu_p) );

    k_sim(t) = k_ss + dr(8,1)*(k_sim(t-1) - k_ss) + dr(8,2)*(n_sim(t-1)     -     n_ss) + dr(8,3)*(l_sim(t-1)-l_ss) ...
                    + dr(8,4)*(rb_sim(t-1)-rb_ss) + dr(8,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(8,6)*(a_sim(t-1)-a_ss);
    n_sim(t) = n_ss + dr(9,1)*(k_sim(t-1) - k_ss) + dr(9,2)*(n_sim(t-1)     -     n_ss) + dr(9,3)*(l_sim(t-1)-l_ss) ...
                    + dr(9,4)*(rb_sim(t-1)-rb_ss) + dr(9,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(9,6)*(a_sim(t-1)-a_ss);
    l_sim(t) = l_ss + dr(10,1)*(k_sim(t-1) - k_ss) + dr(10,2)*(n_sim(t-1)     -     n_ss) + dr(10,3)*(l_sim(t-1)-l_ss) ...
                    + dr(10,4)*(rb_sim(t-1)-rb_ss) + dr(10,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(10,6)*(a_sim(t-1)-a_ss);
    rb_sim(t) = rb_ss + dr(11,1)*(k_sim(t-1) - k_ss) + dr(11,2)*(n_sim(t-1)     -     n_ss) + dr(11,3)*(l_sim(t-1)-l_ss) ...
                      + dr(11,4)*(rb_sim(t-1)-rb_ss) + dr(11,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(11,6)*(a_sim(t-1)-a_ss);
    rkhats_sim(t) = rkhats_ss + dr(12,1)*(k_sim(t-1) - k_ss) + dr(12,2)*(n_sim(t-1)     -     n_ss) + dr(12,3)*(l_sim(t-1)-l_ss) ...
                              + dr(12,4)*(rb_sim(t-1)-rb_ss) + dr(12,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(12,6)*(a_sim(t-1)-a_ss);
    xi_sim(t) = xi_ss + dr(3,1)*(k_sim(t-1) - k_ss) + dr(3,2)*(n_sim(t-1)     -     n_ss) + dr(3,3)*(l_sim(t-1)-l_ss) ...
                      + dr(3,4)*(rb_sim(t-1)-rb_ss) + dr(3,5)*(rkhats_sim(t-1)-rkhats_ss) + dr(3,6)*(a_sim(t-1)-a_ss);
    a_sim(t) = rhoa_p * a_sim(t-1) + shocks(t)*siga_p;

    y_sim(t)   = exp(a_sim(t)) * k_sim(t-1)^alpha_p * h_sim(t)^(1-alpha_p);
    pib_sim(t) = (alpha_p*y_sim(t)/k_sim(t-1)+1-delta_p)*k_sim(t-1) - rb_sim(t-1)*(1-1/l_sim(t-1))*k_sim(t-1) - n_sim(t-1);
    
    if pib_sim(t) > 0
        n_sim_pib(t) = n_sim(t);
%        n_sim(t) = chi1_p*( chi0_p*pib_sim(t) + n_sim(t-1) ) + n0_p;
    else
        minus_profit = minus_profit + 1;
        n_sim_pib(t) = chi1_p*( pib_sim(t) + n_sim(t-1) ) + n0_p;
%        n_sim(t) = chi1_p*(        pib_sim(t) + n_sim(t-1) ) + n0_p;
    end
    
end

minus_profit

return

%% regression to get initial guess for PEA

drop_t = 500;
for t=drop_t:T
    log_rhsee(t-drop_t+1)    = log(lambda_sim(t));
    log_rb(t-drop_t+1)       = log(rb_sim(t));
    log_state_k(t-drop_t+1)  = log(k_sim(t-1));
    log_state_rb(t-drop_t+1) = log(rb_sim(t-1));
    log_state_n(t-drop_t+1)  = log(n_sim(t-1));
    log_xi(t-drop_t+1)       = log(xi_sim(t));
end

state_matrix = horzcat(log_state_k',log_state_rb');
state_matrix = horzcat(state_matrix,log_state_n');
state_matrix = horzcat(state_matrix,log_xi');

reg_rhsee = regstats(log_rhsee,state_matrix);
reg_rb    = regstats(log_rb,state_matrix);

eta_rhsee = reg_rhsee.beta;
eta_rb    = reg_rb.beta;

save('eta_rhsee.txt','eta_rhsee','-ascii','-double');
save('eta_rb.txt','eta_rb','-ascii','-double');

save model_results
